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We study the evolution of magnetized and rigidly rotating neutron stars within a fully general 
relativistic implementation of ideal magnetohydrodynamics with no assumed symmetries in three 
spatial dimensions. The stars are modeled as rotating, magnetized polytropic stars and we examine 
diverse scenarios to study their dynamics and stability properties. In particular we concentrate 
on the stability of the stars and possible critical behavior. In addition to their intrinsic physical 
significance, we use these evolutions as further tests of our implementation which incorporates new 
developments to handle magnetized systems. 
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I. INTRODUCTION 

Neutron stars play a key role in some of the most in- 
teresting astrophysical events observed, from supernova 
remnants and pulsars to a less certain role in long gamma 
ray bursts (GRBs). As such they have attracted signif- 
icant research into their dynamics (for a recent review 
see [l[). In this paper, we revisit the dynamics of rotat- 
ing, and possibly magnetized, neutron stars modeled as 
polytropic stars within a fully nonlinear, general relativis- 
tic model of ideal magnetohydrodynamics (GR+MHD). 
We study the stability properties of these stars and high- 
light possible critical behavior exhibited by the system. 
Furthermore, our studies serve to demonstrate the effec- 
tiveness of our code and certain new developments dis- 
cussed here. This code has been applied to other astro- 
physical problems [l-Q] . 

In recent years, a number of fully relativistic evolutions 
(as opposed to those using Newtonian gravity or other ap- 
proximations to general relativity) of rotating polytropcs 
have appeared studying gravitational wave production 
and black hole formation in astrophysically relevant sys- 
tems. To this end, studies beyond spherical symmetry 
are required which are computationally more demand- 
ing. For instance, rotating stars require moving beyond 
spherical symmetry and many interesting axisymmetric 
scenarios can be addressed using 2D implementations al- 
lowing for excellent resolution without requiring major 
resources (e.g. [§4l0|). In contrast, studies of the most 
general flows and instabilities require 3D simulations and 
are the most expensive realistically accessible scenarios 
being currently considered. Future efforts, including ra- 
diation transport mechanisms, will move beyond 3D sce- 
narios and require efficient use of pctaflop (and beyond) 
resources (e.g. 11|). 

Another important distinction in these studies is 
whether the modeled stars are rigidly rotating or allow 
for differential rotation. Rigidly rotating stars support 
more mass than non-rotating ones (so-called TOV stars), 



but generally do not support the largest masses achieved 
with differential rotation nor do they demonstrate some 
of the more interesting instabilities such as the bar mode 
instability |l2| . Instead, uniformly rotating stars tend to 
demonstrate one of two behaviors; stability or instability 
to collapse to a black hole. Previous studies suggest that 
significant disks do not form from such collapse [l3T - fl5l ]. 
Recall that stars rotating differentially are expected to 
settle into rigidly rotating configurations on short time 
scales, and hence a normal neutron star observed today 
is generally expected to be rigidly rotating. 

Achieving ever more realism, models of neutron stars 
have also begun to consider stars with magnetic field [t| 
Magnetic fields provide, in particular, an effec- 
tive pressure which generally supports greater mass [lj| 
as well as an efficient way to transport angular momen- 
tum. In differentially rotating stars, even small mag- 
netic fields can be amplified by the magnetorotational 
instability. Generally these studies begin with a non- 
magnetized neutron star to which a small, seed magnetic 
field is added. However, in (20[, fully consistent, magne- 
tized stars are evolved, although only nonrotating results 
are presented. Another approach is to study the modes 
through a perturbation approach and analyze the 

growth of these modes with respect to a given stationary 
solution. The dynamical behavior of magnetized stars is 
important also for their role in explaining strong electro- 
magnetic emissions. Indeed, isolated stars with strong 
magnetic fields, so called magnetars, are suspected to be 
the engines powering anomalous X-ray pulsars (AXPs) 
and soft gamma ray repeaters (SGRs) [23j. At least 10% 
of all neutron stars [24j are born as magnetars. In the 
context of single stars, it is thus interesting to exam- 
ine if strong magnetic fields may deform stars away from 
axisymmetry making them strong producers of gravita- 
tional waves 
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Here we present results with uniformly rotating neu- 
tron stars which possess a fully consistent magnetic 
dipole moment. That is to say, the initial data used here 
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represents a stationary state of the full GR+MHD equa- 
tions. These evolutions are computed with a general rel- 
ativity code employing the generalized harmonic scheme 
(allowing for black hole excision). Further, no symme- 
try assumptions are made. High resolutions are achieved 
using a distributed adaptive mesh refinement infrastruc- 
ture. These evolutions demonstrate that our code can 
evolve a stable rotating star for many periods accurately 
Similarly unstable stars evolve to black holes with no 
evidence of any significant disk forming. Finally, we give 
evidence that unstable, rotating, magnetized stars rep- 
resent minimally unstable solutions which could serve as 
Type I critical solutions. 

In Section |n] we provide details about the formulation 
of the equations. In Section ITTTI we discuss aspects of our 
numerical implementation, and describe the diagnostic 
quantities evaluated in Section HVl In Section [V] we dis- 
cuss the initial data we use. We present our results in 
Section [VTT and conclude in Section IvTlI 

II. FORMULATION AND EQUATIONS OF 
MOTION 

Neutron stars can be modeled by rclativistic fluids 
(possibly with the inclusion of magnetic fields) under the 
action of strong gravitational fields (2(| . These systems 
are governed both by the Einstein equations for the ge- 
ometry and by the relativistic equations of magnetohy- 
drodynamics for the matter. We write both systems as 
first order hyperbolic equations. This form of the equa- 
tions is convenient in order to take advantage of several 
rigorous numerical techniques devised for such systems 
to ensure, at the linear level, stability of the implementa- 
tion. More information regarding the motivation for this 
approach can be found in [2|, H^l HH ■ By way of notation, 
we use letters from the beginning of the alphabet (a, b, 
c) for spacetime indices, while letters from the middle 
of the alphabet (i, j, k) range over spatial components. 
We adopt geometric units where c = G = 1. However, 
as discussed in Section [ill! when appropriate, we rescale 
the value of G to achieve improved accuracy in the con- 
servative to primitive variable conversion stage. 

A. The Einstein equations 

The Einstein equations can be written as a system of 
ten nonlinear partial differential equations for the space- 
time metric g a b- The harmonic formulation of the Ein- 
stein equations exploits the fact that the coordinates x a 
can be chosen satisfying the generalized harmonic condi- 
tion mm 

V c V c z a = -r„ = H a (t,x l ), (1) 

where T a = g bc T abc are the contracted Christoffel sym- 
bols. The arbitrary source functions H a (t,x l ) deter- 
mine the coordinate freedom of Einstein equations. The 



original harmonic coordinates correspond to the case 
H a (t,x l ) = 0, which is the choice here. The Einstein 
equations can be expressed in their generalized harmonic 
form [29j j . in particular 

g cd d cd g a b + d a H b + d b H a = -16 n (r ab - j g ab ^j 

+2 T cab H c + 2 g cd g ef (d e9ac d f g bd - Y ace T bd ^j . (2) 

The matter is coupled to the geometry by means of the 
stress energy tensor T ab and its trace T = g ab T ab , which 
will be dictated by the particular model of magnetized 
fluid under consideration, detailed in the next subsection. 

The spacetime can be foliated into hypersurfaces of 
constant coordinate time x Q = t = const. On these 
spacelike hypersurfaces, one defines a spatial 3-metric 
hij = gij. A vector normal to the hypersurfaces is given 
by n a = — V a £/||V a f||, and coordinates defined on neigh- 
boring hypersurfaces can be related through the lapse, a, 
and shift vector, /3\ With these definitions, the space- 
time differential element can then be written as 

d.s 2 = g ab dx a dx b 

= -a 2 dt 2 + h i:i (dx l + f3 l dt) (dx 3 + f3 3 dt) . (3) 

Indices on spacetime quantities are raised and lowered 
with the 4-metric, g ab , and its inverse, while the 3-metric 
hij and its inverse are used to raise and lower indices on 
spatial quantities. 

We adopt a first order reduction of the second order 
differential equations represented in Eqs. <j2j) . This re- 
duction can be achieved by introducing new independent 
variables related to the time and space derivatives of the 
fields 

Qab = -n c d c g ab , Diab = d t g ab ■ (4) 

Within these definitions we can write our evolution equa- 
tions in our GH formalism in the following way [3lj 

d t g a b = P k D kab - a Q a b, (5) 
dtQab = /3 k d k Q ab - ah l3 d l D jab 

- a d a H b - a d b H a + 2 a T cab H c 

+ 2 ag cd (h 13 D lca D jdb - Q ca Qdb - g ef F ace r bd f) 

- ^n c n d Q cd Qa b - a h rj D iab Qj c n c 

- 8ir a(2T ab - g ab T) 

- 2cr a [n a Z b + n b Z a - g ab n c Z c ] 

+ <Ti[3 l (D ulb -d igab ), (6) 

d t D iab = (3 k d k D iab - a d l Q ab 

+ -n c n d D icd Q ab + a h 3k n c D lJC D kab 

- £7i a {Diab ~ dtgab)- (7) 

This GH formulation includes a number of constraints 
that must be satisfied for consistency. On one hand, 
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there are two sets of first order constraints, obtained from 
Eqs. © and defined as 



Ciab 



dig a b — D ia b 
diDj a b 







djDiab 



(8) 



which were introduced when performing the reduction 
to first order 0, Hl[. On the other hand, there are 
the Hamiltonian and momentum constraints, that in the 
Generalized Harmonic formulation show up in terms of a 
four- vector Z a , which is defined as 



2Z a = -T a - H a {t,x l ) 



(9) 



It can be shown that the Hamiltonian and momentum 
constraints are satisfied if Z a = dtZ a — [32|. In order 



to dynamically control the violation of the constraints, 
we have included certain terms proportional to these con- 
straints © These additional terms depend on free 
parameters ao and o\ , allowing one to dynamically damp 
constraint — including the Hamiltonian, momentum, and 
first-order constraints ((8]) — violating modes on a time 
scale proportional to — Oi (|3ll[33l|). 

We evolve the gravitational field equations shown in 
Eqs. ([5][7]) . These equations rely on the computation of 
the 4-dimcnsional Christoffel symbols from the metric g a b 
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(10) 



While we evolve the Di a b functions, the set D$ a b are not 
evolved, but are calculated from evolved quantities as 



Doab 



-aQab + P D kab . 



(11) 



This description suffices to explain the gravitational 
evolution, and the following section describes the evolu- 
tion of the matter. However, we note here that the MHD 
equations are written in the standard 3+1 decomposi- 
tion of spacetime and thus require the spatial metric hij , 
the lapse a, shift f3 l , and ADM extrinsic curvature, Kij. 
These quantities can be written in terms of our evolved 
fields using 



Ki 



= 9ij , a = ^T/g^ , P = 
= -Qij ■■ + — (-D(ij)o - P k D(ij)k 



(12) 



Conversely, the Hamiltonian and momentum constraints 
are usually written in terms of spatial derivatives of the 
metric D^ij and the extrinsic curvature Kij. In fact, we 
use these 3+1 quantities (and similar expressions for their 
derivatives) to calculate the residuals of the Hamiltonian 
and momentum constraints expressed in their standard 
form. 



B. MHD equations 



The stress-energy tensor for the perfect fluid in the 
presence of a Maxwell field is given by 



Tab = [p (l + e) + P}u a Ub + Pg a b 
+ F a c Fb c — -g a bF cd F c d . 



(13) 



The fluid is described by rest mass density p , the specific 
internal energy density e, the isotropic pressure P and 
the four velocity of the fluid u a . With these quantities 
we can construct the enthalpy 



h e = Po + Po( + P, 



(14) 



and construct the standard spatial coordinate velocity of 
the fluid v l as 



(15) 



where W is the Lorentz factor between the fluid frame 
and the fiducial ADM observers. 

The Maxwell tensor F a b can be written as 

F ab = n a E b _ n b E a + fbed ^ ^ ^ 

where E a and B a are the electric and magnetic fields 
measured by a "normal" observer n a . Consequently, both 
fields are purely spatial, i.e., E a n a = B a n a = 0. 

The evolution of the magnetized fluid is described by 
different sets of conservation laws. The magnetic field, in 
the ideal MHD limit, follows the Maxwell equation 



(17) 



where *F ab = e abcd F cd /2 is the dual of the Maxwell ten- 
sor and we have introduced a real scalar field ^ to control 
the divergence constraint. This technique is known as di- 
vergence cleaning [36| and allows for a convenient way to 
control the constraint violation by inducing a damped 
wave equation for the scalar field 'J. The other Maxwell 
equation, in the ideal MHD limit, only gives the definition 
for the current density, since the electric field is given in 
terms of the velocity of the fluid and the magnetic field, 
that is, 



Bk = -e l]k v k B k . 
Conservation of the stress energy tensor in Eq. ([13 
V a T Qb = 0, 



(18) 



(19) 



provides 4 evolution equations for the fluid variables, 
namely the velocity and the internal energy. Conserva- 
tion of the baryon number 



Va(PoU a ) = 



(20) 



We now briefly introduce the perfect fluid equations. 
Additional information can be found in our previous 
work [27], [28| as well as in topical review articles [13, [H| . 



leads to the evolution equation of the rest mass density 
p . Closure of the equations is achieved by introducing 
an equation of state (EOS) relating the pressure with 
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the other thermodynamical quantities of the fluid, P = 
P(Po,e). 

High resolution shock capturing schemes (HRSC) are 
robust numerical methods for compressible fluid dy- 
namics. These methods, based on Godunov's seminal 
work [37j . are fundamentally based on expressing the 
fluid equations as integral conservation laws. To this end, 
we introduce conservative variables q = (D, S^r, £T) T , 
where 



D = W Po , 
Si = (h e W 2 
t = h e W 2 
1 

~2 



I- B 2 )vi — {B^Vj)Bi, 
B 2 -P 

W 2 



D. 



(21) 
(22) 

(23) 



and B 1 is both a primitive and conservative variable. 
In an asymptotically flat spacetime these quantities are 
conserved, and are related to the total energy, momen- 
tum, and, in the non-relativistic limit, the kinetic en- 
ergy, respectively. The quantities w = (p ,v l , P, B l ) T 
are called the primitive variables in contrast to the con- 
servative variables. The fluid state can be specified using 
either set of variables, and both sets are required to write 
the MHD evolution equations. Anticipating the form of 
these equations, we also introduce the densitized con- 
served variables 



D = VhD, S l = VhS i 



B' 



VhB\ 



(24) 

where h = det(ft.y). The fluid equations can now be 
written in balance law form 



d t q + d k f k (q) = s(q) 



(25) 



where f k are flux functions, and s are source terms. The 
fluid equations in this form are 



d t D + di 



dtSj + 8, 



aD\v l - — 
a 



= 0, 



(26) 



Sj (v l -—) +VhPh l 3 



« 3 r' 



a S l 



[s iV k + Vh~j 

-dja(f + D) 
-(a(B t W~ 2 



Phf + s a d^ a 



ViVj&)d k B k , (27) 



v l D 



1 X, 



K ii S i v j +VhKP-- S a d a a 



-Cav J B 1 d k B k 



d t B b + 8 t 



B b I v 



a 



dt$> = -c r a^ 



(28) 



(29) 



Ch-^diB 1 + (F .(30) 



where c r = k and we have allowed for different speeds 
than light by introducing the parameter cu- Except in 
the tests, we will use the original prescription (fT7|) with 

Ch = 1. Here 3 T l j k is the Christoffel symbol associated 
with the 3- metric /ly, and K is the trace of the extrin- 
sic curvature, K = K l i. Notice that the aforementioned 
system is an extended version of the one employed in our 
earlier works [H, H3, [28| Here we have added additional 
terms toggled by the parameter £ which allow for consid- 
ering an extension of the "eight-wave" formulation which, 
in the absence of the cleaning field VP, ensures the strong 
hyperbolicity of the system [36[ -at the cost of introduc- 
ing derivative terms in the sources-. Furthermore, by 
setting Q = 1, the propagation speeds of two constraint 
violating modes become non-vanishing and hence viola- 
tions are dragged along by the fluid's velocity. This is 
numerically convenient as possible violations will propa- 
gate off the grid. 

Finally, we close the system of fluid equations with an 
equation of state (EOS). We choose the ideal fluid EOS 



P=(T-l)p e, 



(31) 



where T is the constant adiabatic exponent. Nuclear mat- 
ter in neutron stars is relatively stiff, and we set T = 2 
in this work. When the fluid flow is adiabatic, this EOS 
reduces to the well known polytropic EOS 



P 



np 



(32) 



where k is a dimensional constant. We use the polytropic 
EOS only for setting initial data. 

The set of fluid equations, Eqs. (|26H29[) . are used to 
evolve the conservative variables. However, these equa- 
tions also contain the primitive variables which neces- 
sitates a step in the evolution scheme which solves for 
the primitive variables in terms of the conservative ones. 
Given the primitive variables, the conservative variables 
are easily calculated from the algebraic expressions in 
Eqs. (j2~4"]) . Calculating the primitive variables from the 
conservative variables, however, is more delicate, as it re- 
quires the solution of a transcendental equation. To this 
end we first define the quantity x = h e W 2 . We then write 
S 2 = S l St in terms of x and solve for W 2 , obtaining 



W 



x 2 {x + B 2 ) 



2 i2 



! (x + B 2 ) 2 - x 2 S 2 - (2x + B 2 )(S t B 1 ) 2 



(33) 



Using the definition of r, we define a function that is 
identically zero 



f( x ) = x -P-l 



B^_ 

W 2 



-B 2 -t-D = 0. (34) 



If the enthalpy can be expressed as a simple function of 
the pressure, as can be done for the ideal fluid and a 
generalized EOS with a cold nuclear component, then we 
can express the pressure as a function x. For the ideal 
fluid EOS used here, the enthalpy equation 



x 

W 2 



Po + Poe + P 



(35) 
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FIG. 1: Demonstration of the effectiveness of divergence 
cleaning at handling deviations from a divergenceless mag- 
netic field. A calculation of the divergence of the magnetic 

field, x 2 |V • -Bj , along the i-axis for three times is shown for 

each of three cases: (top) No divergence cleaning; (middle) 
Cleaning with Ch = 0.1 and c r = 0.01; and (bottom) Cleaning 
with Ch = 1.0 and c r = 0.1. These runs were otherwise iden- 
tical for a magnetized, rotating star of coordinate equatorial 
radius of 10 with spin along the z-axis with three levels of re- 
finement (at ±50, ±25, and ±12.5) and a perturbation to the 
magnetic field to introduce an explicit deviation from diver- 
genceless. Concentrating on the MHD equations, the metric 
terms were frozen at their initial values, e.g. the Cowling ap- 
proximation. The bottom row shows clear wavelike behavior 
as it "cleans" the divergence. 



can be solved for P by substituting in po 
the EOS to obtain 



P 



T - 1 



x 



D 
W 



D/W and 



(36) 



Combining these equations, / is a function of a single 
unknown x. This equation can be solved numerically 
using the Newton-Raphson method. It is useful to note 
that x has a minimum physical value, which is found by 
requiring in Eq. (|33|) that W 2 > 1. 



III. IMPLEMENTATION DETAILS 

The code is constructed within the HAD computa- 
tional infrastructure which provides distributed adaptive 
mesh refinement (AMR). The AMR follows in the style 



of Berger & Oliger J38(, but uses the tapering condi- 
tion for AMR boundaries instead of temporal interpo- 
lation [HI . The combined set of geometric equations and 
fluid equations, Eqs. ([5][7]) and Eqs. (|26H30p respectively, 
is discretized using the method of lines. The geometric 
equations are discretized using operators that satisfy a 
summation by parts property [40 . |41[ . The fluid equa- 
tions are discretized using the HLLE method 42 1. The 
semi-discrete equations are solved using a third order ac- 
curate, total variation diminishing (TVD) Runge-Kutta 
solver (H. 

The fluid equations diverge as the density goes to zero, 
and, as is standard practice, we disallow any true vacuum 
by setting such regions to a floor or atmosphere value. 
The floor is applied after each fluid update as 



L) 



i(A S), 



min(r, 5), 



(37) 



where S is chosen to be many orders of magnitude smaller 
than the maximum densities and pressures in the initial 
data. The comparison of otherwise identical runs but 
with different floor values suggest that the use of an at- 
mosphere generally does not affect accuracy. 

We have, however, found certain issues with precision 
occurring within the primitive solver. Typical maximum 
values of the density are about 10 -2 in geometric units, 
with a floor value of S = 10~ 8 . We have found it use- 
ful therefore to scale Newton's constant G such that 
the fluid densities and pressures are close to order unity. 
Thus, rather than using the typical choice of G = 1, we 
might use G = 1/1000. As G affects only the coupling 
of the fluid to the geometry, the evolution of the geomet- 
ric equations is not affected by this scaling. Empirically, 
we find that scaling G allows the primitive variables to 
be more easily recovered in low density regions. This im- 
provement appears to be related to finite precision effects 
in the primitive solver. The scaling decreases the effec- 
tive floor value while avoiding the problems associated 
with having a true vacuum. 

The Maxwell equations require that the magnetic field 
be divergenceless. This is the so-called "no monopolc" 
constraint. A variety of schemes exist with the goal 
of controlling the growth of the divergence. We choose 
a strategy that ensures flexibility and robustness when 
dealing with multiple grid structure (as in AMR) and al- 
lows, in principle, for a clean boundary treatment [44| . 
To this end we have implemented hyperbolic divergence 
cleaning as described in [36[ (also see [45j]). We thus 
introduce a scalar ^(x, y, z, t) which is sourced by the 
negative of the divergence of the magnetic field as shown 
in Eq. (|3"0"|) . As described in [36|, this scheme implies 
that the divergence obeys a damped wave equation so 
that constraint violations propagate off the grid and their 
value is reduced. 

For initial data generated with Magstar, the diver- 
gence is around machine precision, and so to test the 
implementation of divergence cleaning, we introduce a 
perturbation to the magnetic field in order to produce a 
significant amount of divergence to the magnetic field. In 
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particular, this perturbation takes the form of a spher- 
ical, Gaussian shell of radius r , width 5, and ampli- 
tude A added to each component of B. We expect the 
divergence cleaning to propagate this perturbation as a 
damped wave, and we therefore plot the scaled quantity 

x 2 V • B in Fig. [TJ As can be seen by comparing the 

cases of no cleaning {top) and with cleaning (bottom), the 
divergence propagates with damping through the refine- 
ment boundaries across the grid. 

As typical with codes dealing with linearly degenerate 
hyperbolic systems, like those those employed in numeri- 
cal relativity, a dissipation operator is applied to the met- 
ric variables. This operation uses a high-order derivative 
to serve as a low-pass filter and does not affect the accu- 
racy of the simulation. We have found it useful for keep- 
ing things smooth. Although this operation is not applied 
to the fluid variables, we have found it quite important 
in keeping the magnetic field components smooth. The 
magnetic field evolution is coupled tightly with that of 
the velocity, and any nonsmoothness which appears in 
the velocity can easily affect the smoothness of the mag- 
netic field. The addition of dissipation to the magnetic 
field and divergence cleaning field helps control the be- 
havior of the magnetic field. 

At the boundaries of the domain, simple outflow 
boundary conditions are applied to the fluid variables. 
This is accomplished by copying the values of the conser- 
vative variables near the boundary outward. Most of the 
gravitational variables are treated using Sommcrfcld-likc 
boundary conditions of the form p6| 



d t + d r + - ) 







(38) 



where rj a b is just the Minkowski metric. The rest, which 
are not so crucial, are set either by maximally dissipa- 
tive or constraint preserving boundary conditions [3lj . 

A relevant issue when considering the collapse of stars 
is the formation of a black hole. To deal with such situa- 
tions, we adopt black hole excision where we dynamically 
monitor for the appearance of trapped surfaces (which lie 
inside an event horizon if cosmic censorship holds) and 
excise cubical region(s) from the computational domain. 
As discussed in [47I . |48| . this excision introduces inner 
boundaries which are of "outflow" type and so no bound- 
ary condition is required there. However, for a more ro- 
bust handling of the fluid, we also allow for a modification 
of the fluid equations inside the trapped surface Q. The 
MHD equations are written in balance law form 



U + F(U)' = S, 



(39) 



which we modify to include a damping term near the 
black hole 



U + F(U)' = S- f(r)(Ax) p (U - U ). 



(40) 



Here the function f(r) decreases smoothly with r. from a 
given value at the excision region to zero at the outermost 



trapped surface (OTS) found, and is zero for r > ^otSj 
so that the exterior of the BH is causally disconnected 
from the effect of this extra term. Uq is set to zero or 
to the value of the atmosphere if the corresponding field 
has one. The coefficient (Ax) p ensures that the damping 
term converges to zero as the grid spacing Ax is reduced. 
As long as one chooses p greater than or equal to the order 
of convergence of the code, this term will not modify the 
convergence rate. We typically adopt a value of p = 4. 

Finally, gravitational radiation is calculated via the 
evaluation of the Newman-Penrose scalar ^4 which is 
computed by contracting the Weyl tensor, C a tcd, with a 
suitably defined null tetrad {£, n, to, to} 



*4 = C abcd n a m b n c m a 



(41) 



extracted at spherical surfaces £j located in the wave- 
zone, far from the sources. We also consider possible 
corrections required to deal with gauge ambiguities, as 
discussed in [49j]. We refer the reader to that paper for 
details on the adopted tetrad and required corrections. 



IV. DIAGNOSTIC QUANTITIES 

The initial configuration for the stars is axisymmetric, 
and we therefore want to be able to measure any change 
to this structure. For this purpose, we monitor certain 
distortion parameters [5(| defined as 



jyy 



n+ 
'// 
v 



Jxx _|_ Jyy 

2I xy 

jxx _|_ jyy 



in terms of the moment of inertia tensor 



Dx j x & 



x. 



(42) 
(43) 
(44) 

(45) 



These parameters are computed with respect to the con- 
servative variable D. 

It is also standard practice to display the maximum 
of the (primitive) rest mass density, and we compute the 
fractional change in time as 



Ap. 



maxj£oft)| - max I po (0)| 
max \po(0)\ 



(46) 



Generally, for stable stars one sees this quantity oscillate 
from inherent numerical perturbations about a stable so- 
lution. For unstable solutions, one expects this quantity 
to change in a significant way. Similarly, we compute the 
relative change in baryon mass as a function of time as 



-^baryon(^) -^baryon (0) 



Mi 



baryo 



n(0) 



(47) 
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where 

M baryon = J DdV. (48) 

This mass is related to the expected baryon number and 
should be strictly conserved as long as mass is not leav- 
ing the computational domain. Many HRSC schemes 
explicitly conserve this quantity, but here it is not a pri- 
ori conserved. Our use of an atmosphere entails adding 
a small amount of mass in regions that would otherwise 
become evacuated. Another reason is that we are using 
a finite difference based AMR which does not accom- 
modate as readily a strictly conservative treatment as a 
finite volume method would (HT| ■ Finally, the presence 
of source terms in the evolution equations for the other 
conservative variables also breaks perfect conservation of 
these other quantities. 

We also compute the angular momentum of the fluid. 
Since our stars rotate about the z-axis, we need only 
compute a single such quantity 

J z = J {xn j T ]y - yn j T JX ) Vh d 3 x. (49) 

The use of Cartesian coordinates is typically avoided 
when one expects angular momentum to be conserved 
since cylindrical coordinates are better adapted to the 
respective Killing vector. However, because this project 
is part of a more general effort to model a variety of 
systems with different natural coordinates and no sym- 
metries, we simply monitor the extent to which this is 
conserved by computing a fractional change as 

Finally, we also monitor the extent to which our nu- 
merical solution satisfies the Einstein equations as ex- 
pressed in Eq. (fTJ. That is, we compute the so-called 
residual of these equations, expressed as the four-vector 
Z a as in Eq. ([9]). In particular, we monitor the norm of 
this vector 

\Z\ = \\Z a \\2. (51) 

Analytical solutions to Einstein equations must have a 
vanishing residual and thus we monitor this residual for 
solutions obtained at various resolutions to check for con- 
vergence to a consistent solution. 

V. INITIAL DATA 

We use initial data generated with the program 
Magstar, part of the Lorene software package. These 
solutions are described in [lj| and are rigidly rotating, 
magnetized neutron stars. They are generated as fully 
consistent solutions of Einstein's equations as opposed 
to taking nonmagnetized stars and adding a seed mag- 
netic field without re-solving the constraints. The initial 




0.2 0.4 0.6 0.8 1 

H 

c 

FIG. 2: Diagram of solution space of stars generated with 
Magstar (using its units). The total gravitational mass is 
plotted versus the central enthalpy. The upper curve (pur- 
ple, open triangles) represents stars with no magnetic field 
at the mass shredding limit. Very close to this are shown 
(green, open circles) stars rotating at the same frequencies but 
with a radial magnetic field at the pole of 1000 GT = 10 16 G. 
The masses of these stars are barely changed with respect 
to their nonmagnetic counterparts. The lower curve (blue, 
open squares) represents the static limit with no rotation. A 
sequence of constant baryon mass stars (Mb = 3.6M ao i) is 
shown (red stars). This sequence is also shown in the inset 
along with their angular momentum (yellow open circles) and 
frequency of rotation (green stars). The location of the initial 
data used in the initial data convergence test of Fig. [3] and 
the long term, stability test of Fig. [4] is also shown (green, 
solid square) although that solution strictly does not belong 
here because it has non- vanishing magnetic field. 



magnetic field is dipolar and aligned with the rotation 
axis, produced by a current function f{A<p) = constant 
where is the toroidal component of the electromag- 
netic potential vector. 

To get a handle on the solutions generated, we first 
turn off the magnetic field and compute the "usual" two- 
parameter solution space as described in [52| ■ Using the 
terminology of (l9l |. we compute solutions based on the 
two parameters of central enthalpy and frequency. Ex- 
amination of Eq. (13) of [lH shows that the log-central 
enthalpy H c is related to the entropy used here (h e as 
defined in Eq. ([13])) by 

H c = In (h e ) + C (52) 

where C is constructed by physical constants and h e is 
evaluated at the center of the star. 
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FIG. 3: The residuals of two constraint equations, namely 
the Hamiltonian constraint and the y-component of the mo- 
mentum constraint. A slice along the z-axis of the com- 
puted residuals for various unigrid resolutions. Top: That 
the Hamiltonian constraint residual converges to zero with in- 
creasing resolution is taken as evidence that the initial data is 
properly constructed and read into the code. Bottom: That 
the momentum constraint residual improves with the first in- 
crease in resolution is similar evidence. However, the next 
increase in resolution fails to bring down this residual. We 
suspect that this remaining error is associated not with trun- 
cation error but instead with inherent errors in the numeri- 
cal transformation and interpolation from Lorene's spherical 
basis to our Cartesian one. Note the spikes that arise at the 
stellar boundaries due to discontinuities in the fluid variables, 
these are not expected to converge to zero. 



quences it has been shown that there is a stability change 
at the minimum in the angular momentum [52j | . Looking 
at the inset of Fig. [51 one therefore expects the solutions 
on the right to be unstable while those on the left should 
remain stable. 

Note that the addition of a non-vanishing magnetic 
field adds another dimension to this diagram although 
the effect of the magnetic field on the initial data is not 
so dramatic. Certainly there are many ways in which to 
"add" a magnetic field to a stellar solution. In [l9| a non- 
vanishing current is assumed and a solution is obtained 
with the same baryon mass but now with nonvanishing 
magnetic field. In particular, solutions at the mass shed- 
ding limit with no magnetic field arc shown in Fig. [5] by 
open triangles. Almost indistinguishable from these solu- 
tions are the those magnetized such that the radial mag- 
netic field at the pole of the star is 1000 gigatesla (GT) 
or 10 16 G. shown with open circles. 

We consider different perturbations. To perturb the 
pressure, we decrease it according to real parameters A p 
and m in terms of the unperturbed pressure po depending 
on the azimuthal angle tp 



V 



Po[l 



A p sin 



(imp)] . 



(53) 



We also consider perturbations to the rest mass density 
decreasing it everywhere by a fraction A p 



p = p [l- A p ] 



(54) 



We verify the initial data in our evolution code in a 
number of ways. We examine convergence of the data to a 
unique solution which solves the constraints. In Fig. [3] we 
show the residual of the Hamiltonian and y-component of 
the momentum constraint for three different resolutions. 
The Hamiltonian constraint is a good measure of the fi- 
delity of the solution to the Einstein equations, and, as is 
clear from the figure, its residual decreases rapidly with 
resolution. Furthermore, we confirm that the divergence 
of the magnetic field is around machine precision. 



These nonmagnetized solutions are diagrammed in a 
plot of total gravitational mass versus central enthalpy 
in Fig. [2j The lower curve shows the static limit for stars 
which are not rotating. The upper curve is the mass 
shedding limit, represented by the largest frequency for 
which Magstar returned a solution. These curves serve 
as the upper and lower bounds on the solution space of 
stationary, unmagnetized stars. It should be noted that 
Magstar generates only rigidly rotating stars, and there- 
fore the evolutions are far from the fast rotating regime 
expected to excite large and growing nonaxisymmetric 
modes. 

We also compute and show a sequence of stars at con- 
stant baryon mass. Such sequences are important be- 
cause real stars are expected to conserve baryon mass as 
they evolve and thus the sequences are expected to ap- 
proximate their evolution. Furthermore, along such se- 



VI. RESULTS 
Stable, Rotating Star Test: 

Before addressing the effects of the magnetic field, we 
verify that the code reproduces the expected behavior as 
described, for example, in [53j]. First, we consider evo- 
lutions of stable, rotating stars and find that the code 
evolves such a star as long as desired while maintain- 
ing a stationary solution. This a demanding test as it 
depends on the balance between gravitational and hy- 
drostatic forces in a rotating configuration with both a 
non-conforming grid and variables not adapted to the 
symmetries of the problem. 

One example of the the behavior of the numerical so- 
lution is shown in Fig. [4] As apparent in the top frame of 
the figure, the fractional change in the maximum of the 
density oscillates with a slow overall increase. This oscil- 
lation is characteristic of quasinormal ringing of the star 
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FIG. 4: Convergence results for the evolution of a stable, 
nonmagnetized, rotating star (the particular initial solution 
is shown in Fig. [2] as a green solid square). Three FMR 
evolutions are shown, each with a multiple of the coarsest 
resolution resolving the equator of the star with: (magenta, 
dot-dashed) 30 points/star, (blue, dotted) 60 points/star, and 
(black, solid) 120 points/star. Also shown is a run with the 
same stellar resolution of (red, long-dashed) 60 points/star 
but with a coarse grid that extends twice as far in all direc- 
tions. The two highest resolution runs were terminated early 
because of the computational cost, not because of any robust- 
ness problems. The insets show the same data but in finer 
detail for the first rotational period. The top frame shows 
the fractional change in maximum density versus rotational 
period. With increasing resolution the solution better ap- 
proximates a stationary solution. The upper middle frame 
shows the fractional change in the baryon mass, Although our 
scheme, using vertex centered AMR, an atmosphere, and full 
general relativity with sources, is not strictly conservative, the 
plot shows that deviations from strict conservation are small 
and decrease with more resolution. The lower middle frame 
shows the fractional change in integrated angular momentum. 
Again, the code demonstrates convergence to conservation. 
The bottom frame shows the maximum of the norm of the 
Z a constraint residuals as a function of time which also con- 
verge. 



excited by inherent numerical error. However, the figure 
shows data for a number of resolutions, and the trend as 
resolution improves is toward a flatter curve. This trend 
is particularly apparent in the inset which shows just the 
first rotational period. This behavior suggests that the 
code is converging to the continuum solution. 

In addition to the runs with varying resolution, Fig. |4] 
shows another evolution with a coarse level extending 



FIG. 5: Collapse to black hole of an unstable, unperturbed, 
magnetized star for multiple of a base resolution. The ini- 
tial data is an unperturbed Magstar solution with central en- 
thalpy H c = 0.8 rotating at a frequency / = 835Hz and with 
polar magnetic field of 1000 GT = 10 16 G. The top frame 
shows the maximum density which increases with time as the 
star collapses. The upper middle frame shows a norm of 
the divergence of the magnetic field. The lower frames show 
the distortion parameters of the D field as functions of time. 



twice as far but with fine levels identical to the medium 
resolution run just discussed. Generally, the results are 
the same as for the non-extended domain, indicating little 
effect from the boundary for the the first half-period. 

Although the fluid scheme used here is not strictly con- 
servative because of (i) AMR boundaries for our vertex 
centered scheme, (ii) our outer boundary treatment, and 
(iii) our use of a fluid floor, Fig. @] shows that the frac- 
tional change in the volume-integrated baryon mass re- 
mains constant to a high degree. Similarly, the integrated 
angular momentum of the fluid converges to conserva- 
tion. 

Finally, the bottom frame of Fig. @] shows the maxi- 
mum value of the constraint violation. These values in- 
crease with time as the numerical error accumulates, but 
higher resolution runs demonstrate less violation though 
it saturates due to the intrinsic error of the initial and 
boundary data. 

Unstable, Rotating, Magnetized Star Test: 

Similarly, we evolve a rotating, magnetized star located 
on the unstable side. These stars, perturbed by inherent 
numerical error, collapse to black holes. We see no evi- 
dence of significant disk formation (as in e.g. (l5l. [l6l[54T- 
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l56j). Notice that even though we do not impose any 
type of symmetry, no asymmetric, unstable modes are 
observed. This behavior is consistent with previous work 
which studied the possible onset of axisymmetric insta- 
bilities [57| . These previous studies found that such in- 
stabilities require high rotation rates characterized by 
T/W > 0.25 in contrast to those studied here for which 
T/W < 0.1. Consequently, we see no evidence for defor- 
mations of the star as would be apparent by monitoring 
the distortion parameters moving away from zero. 

In Fig. [5j we show an example of such an evolution at 
three successively finer resolutions. Because the star is 
collapsing, the maximum density increases dramatically. 
The magnetic field remains essentially poloidal through- 
out the collapse and its maximum magnitude increases 
due to the resulting compression of the field lines. De- 
spite this increase, the magnetic field plays no important 
role in the collapse since its associated pressure is still 
several orders of magnitude smaller than the fluid pres- 
sure. While we do observe an increase in the norm of the 
divergence, the growth is not particularly fast and the 
divergence remains small in absolute terms and relative 
to the magnetic field. 

Fig. [5] displays two snapshots of the density and mag- 
netic field strength of the collapsing star tested in Fig. [5] 
The first one illustrates a stage during the collapse before 
an apparent horizon forms. The second one shows the 
behavior after an apparent horizon is found and its inte- 
rior excised. The apparent horizon appears at t ~ 0.6 P, 
when the maximum of the density is p ma x = 0.271 and 
the minimum of the lapse is a m in = 0.2. 

Perturbations of Unstable Stars: 

Previous work presented in [53[ argues that, in general, 
unstable stars should either collapse to a black hole or ex- 
pand and oscillate about a stable stellar solution. Seeking 
to duplicate this behavior, we perturb an unstable star. 
Indeed, as shown in Fig. [7J we find precisely the behav- 
ior described. If our perturbation increases the pressure 
sufficiently, we find a star which expands and oscillates 
about some other, presumably stable solution. However, 
if we choose a perturbation which barely increases the 
pressure, the star collapses to a black hole. 

However, given the interest in black hole critical phe- 
nomena (see [H, HI]) over the past couple of decades, 
we study in detail the separation between these two be- 
haviors. That is, we continue to adjust A p in Eq. (|53[) 
searching for a value A* above which one finds black hole 
formation, and below which one finds an expanding so- 
lution (the way we have parameterized the pressure per- 
turbation, A* < 0). 

As apparent from Fig. [JJ the more one continues this 
tuning, the longer the unstable star survives. This type of 
tuning is reminiscent of a similar analysis of an unstable, 
irregular static solution [60[. What these results suggest 
is that the unstable solution (i) sits at the threshold of 
black hole formation with (ii) a single unstable mode. 
That small perturbations about the solution send it ci- 




t/P 

FIG. 6: Density po and magnetic field strength on the y — 
plane for the same star as shown in Fig. [5] The density is 
shown with respect to the colormap while the contours denote 
the magnetic field, which is equally spaced from to 5x 10 15 G. 
The top left corresponds to t = 0.5P while the the top right to 
t = 0.9P (the circle shown represents the apparent horizon.) 
The bottom plot illustrates the normalized mass as a function 
of time, along with its rate of change. 




t/P 

FIG. 7: The maximum density p max as a function of time for 
an unstable, perturbed (m = 3 in Eq. (|53[) h non-magnetic, 
rotating star. The (unperturbed) initial star has central en- 
thalpy 0.8 and rotates at the mass shedding limit. Solutions 
resulting from tuning the amplitude of the perturbation A p to 
roughly one part in 10 6 show the two disparate outcomes (as 
discussed in [53|]): collapse to black hole or violent oscillations 
about a stable stellar solution with equivalent mass. With 
more tuning, the star resembles the initial, unstable solution 
for a longer time. 
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ther to collapse or expansion suggest that it sits at the 
threshold. Furthermore that a single parameter is suffi- 
cient to stabilize the solution suggests that there is a sin- 
gle unstable mode. In contrast, sometimes one can tune 
multiple parameters to find a threshold solution [61 1. 

If these suggestions hold up, these would suggest that 
at least some of these unstable solutions might serve 
as Type I critical solutions. Indeed, previous work 62 1 
perturbed stable TOV stars in spherical symmetry 
and found unstable TOV stars at criticality in Type 
I collapse. In that work, they were able to achieve 
phenomenal resolution and tuning. In contrast, while 
they perturbed a self-consistent stable solution and saw 
the tuned evolution driven to the unstable branch of 
solutions, here we begin with the unstable solution and 
perturb around it. Here, we have only been able to tune 
to about one part in a million, being prevented from 
tuning further because successive evolutions stop the 
trend towards longer lived solutions. That such searches 
are prevented from continuing might indicate some 
new phenomena, or, more likely, that boundary effects 
and numerical error begin to spoil the threshold behavior. 

We have looked for a scaling law in the survival times of 
these tuned unstable stars. To the extent that this rough 
tuning is representative of the overall behavior, we find 
that the different solutions appear to scale as expected. 
However, because our searches have terminated so far 
from criticality, we cannot have much confidence in a 
precise scaling relationship. There has been recent work 
in the axisymmetric collision of neutron stars [63l I&H 
which appears to demonstrate the same type of critical 
behavior as observed in [62^ . 

We find the same type of behavior about a magne- 
tized star as shown in Fig. [8j Here we have carried out 
three searches by varying a different parameter all per- 
turbing the same solution. The figure makes apparent the 
same ringing for all three families, although the ampli- 
tude varies across the different tuning families. It seems 
reasonable to take the results of these tunings as further 
evidence that only a single mode is unstable since if there 
were more unstable modes, these different tunings would 
produce solutions more varied from each other. We have 
begun to look at the geometry of the purported unstable 
mode by looking at the difference A = pa(t) — po(0) at 
late times for near-critical evolutions. These calculations 
indicate the the mode is likely axisymmetric with differ- 
ences between A on the x = plane and on the y = 
plane being at about the 10% level. 



0.164 
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FIG. 8: Results of three critical searches for a magnetized, ro- 
tating star with central enthalpy of 0.8 and a radial magnetic 
field at the pole of 1,000GT. One search varied the ampli- 
tude of the pressure perturbation A p with m — 7 (black and 
cyan), another with m = 3 (blue and magenta), and another 
varied the density perturbation A p . In solid line are shown the 
super-critical evolutions which collapse to black holes, while 
sub-critical solutions for the two searches are shown with dot- 
ted lines. The search over pressure with m = 7 achieved tun- 
ing to about one part in 10 3 , that with m = 3, one part in 
10 4 , and that over density achieved about one part in 10 5 . 



is part of the publicly available Lorene package. We 
present several studies which indicate our code reliably 
and robustly evolves astrophysically relevant scenarios 
including magnetic field effects. We study the dynamical 
effect of perturbations of stars on the unstable branch of 
solutions. We find evidence that these unstable solutions 
may play a similar role as the unstable TOV st ars p lay in 
spherically symmetric evolutions as studied in |62j . This 
is significant because it suggests that the addition of an- 
gular momentum, magnetic field, and three dimensions 
do not allow for a multitude of unstable modes. Needless 
to say, the phenomena associated with the threshold of 
gravitational collapse merit further study. 

Beyond the studies considered in this work, further in- 
teresting phenomenology to consider include the impact 
of magnetic fields in the stability of the star, a thor- 
ough analysis of the possible critical phenomena observed 
and differences due to more realistic equations of state. 
The investigation of such scenarios will be presented else- 
where. 



VII. CONCLUSIONS 

We study the evolution of rotating stellar configura- 
tions and examine different phenomenology related to 
stability and magnetic field influence in their dynamical 
behavior. The stars are constructed assuming a poly- 
tropic equation of state using the code Magstar which 



Acknowledgments 

We would like to thank J. Novak for his gracious assis- 
tance with Magstar and Lorene. SLL thanks Matthew 
Choptuik and Scott Noble for valuable discussions con- 
cerning the apparent critical behavior. We also would like 
to thank Matthew Anderson, Eric Hirschmann, Miguel 
Megevand, Patrick Motl, Joel Tohline and Travis Gar- 



12 



ret for comments and suggestions. SLL, DN and CP 
thank the Perimeter Institute for Theoretical Physics 
for hospitality where parts of this work were com- 
pleted. This work was supported by the National Sci- 
ence Foundation under grants PHY-0803629 and PHY- 
0653369 to Louisiana State University, CCF-0832966 and 
PHY-0803615 to Brigham Young University, and CCF- 
0833090, PHY-0803624 to Long Island University, and 
NSERC through a Discovery Grant. Research at Perime- 



ter Institute is supported through Industry Canada and 
by the Province of Ontario through the Ministry of Re- 
search & Innovation. This research was also supported in 
part by the National Science Foundation through Tera- 
Grid resources provided by SDSC under allocation award 
PHY-040027. In addition to TeraGrid resources, we 
have employed clusters belonging to the Louisiana Opti- 
cal Network Initiative (LONI), and clusters at LSU and 
BYU. 



N. Stergioulas, Living Rev. Rel. 6, 3 (2003), gr- 
qc/0302034. [25 
C. Palenzuela, I. Olabarrieta, L. Lehner, and S. Liebling, 
Phys. Rev. D75, 064005 (2007), gr-qc/0612067. 
C. Palenzuela, L. Lehner, and S. L. Liebling, Phys. Rev. [26 
D77, 044036 (2008), 0706.2435. 

M. Anderson et al., Phys. Rev. D77, 024006 (2008), [27 
0708.2720. 

M. Anderson et al., Phys. Rev. Lett. 100, 191101 (2008), [28 
0801.4387. 

C. Palenzuela, M. Anderson, L. Lehner, S. L. Liebling, 

and D. Neilsen, Phys. Rev. Lett. 103, 081101 (2009), [29 
0905.1121. [30 
M. Megevand et al., Phys. Rev. D80, 024012 (2009), [31 
0905.3390. 

M. Shibata and Y.-i. Sekiguchi, Phys. Rev. D69, 084024 
(2004), gr-qc/0402040. [32 
M. D. Duez, Y. T. Liu, S. L. Shapiro, and M. Shibata, 
Phys. Rev. D73, 104015 (2006), astro-ph/0605331. [33 
K. Kiuchi, M. Shibata, and S. Yoshida, Phys. Rev. D78, 
024029 (2008), 0805.2712. [34 

D. Bader, ed., Petascale Computing: Algorithms and Ap- 
plications (Chapman & Hall/CRC, 2007). [35 
B. Zink et al., Phys. Rev. D76, 024019 (2007), astro- [36 
ph/0611601. 

M. Shibata, Astrophys. J. 595, 992 (2003), astro- 
ph/0310020. [37 
M. D. Duez, S. L. Shapiro, and H.-J. Yo, Phys. Rev. [38 
D69, 104016 (2004), gr-qc/0401076. 

L. Baiotti et al., Phys. Rev. D71, 024035 (2005), gr- [39 
qc/0403029. 

B. C. Stephens, S. L. Shapiro, and Y. T. Liu, Phys. Rev. [40 
D77, 044001 (2008), 0802.0200. [41 
M. Shibata, M. D. Duez, Y. T. Liu, S. L. Shapiro, and 
B. C. Stephens, Phys. Rev. Lett. 96, 031102 (2006), 
astro-ph/0511142. [42 
M. Shibata, Y. T. Liu, S. L. Shapiro, and B. C. Stephens, 
Phys. Rev. D74, 104026 (2006), astro-ph/0610840. [43 
M. Bocquet, S. Bonazzola, E. Gourgoulhon, and J. No- 
vak, A&A 301, 757 (1995), arXiv:gr-qc/9503044. [44 
B. Giacomazzo and L. Rezzolla, Class. Quant. Grav. 24, 
S235 (2007), gr-qc/0701109. [45 
H. Sotani, K. D. Kokkotas, N. Stergioulas, and [46 
M. Vavoulidis (2006), astro-ph/0611666. 
H. Sotani, K. D. Kokkotas, and N. Stergioulas, Mon. Not. [47 
Roy. Astron. Soc. 375, 261 (2007), astro-ph/0608626. 
S. Mereghetti, Astron. Astrophys. Rev. 15, 225 (2008), [ 
0804.0250. 

A. M. Beloborodov and C. Thompson, Astrophys. J. 657, 



967 (2007), astro-ph/0602417. 

B. Haskell, L. Samuelsson, K. Glampedakis, and N. An- 
dersson, Mon. Not. Roy. Astron. Soc. 385, 531 (2008), 
0705.1780. 

N. Andersson and G. L. Comer, Living Rev. Rel. 10, 1 
(2005), gr-qc/0605010. 

D. Neilsen, E. W. Hirschmann, and R. S. Millward, Class. 
Quantum Grav. 23, S505 (2006). 

M. Anderson, E. Hirschmann, S. L. Liebling, and 
D. Neilsen, Class. Quant. Grav. 23, 6503 (2006), gr- 
qc/0605102. 

H. Friedrich, Commun. Math. Phys. 100, 525 (1985). 

D. Garfinkle, Phys. Rev. D65, 044029 (2002). 

L. Lindblom, M. A. Scheel, L. E. Kidder, R. Owen, 

and O. Rinne, Class. Quant. Grav. 23, S447 (2006), gr- 

qc/0512093. 

C. Bona, T. Ledvinka, C. Palenzuela, and M. Zacek, 
Phys. Rev. D66, 084013 (2002). 

C. Gundlach, J. Garcia, G. Calabrese, and I. Hinder, 
Class. Quant. Grav. 22, 3767 (2005), gr-qc/0504114. 
J. M. Marti and E. Miiller, Living Rev. Rel. 2, 3 (1999), 
astro-ph/9906333. 

J. A. Font, Living Rev. Rel. 3, 2 (2000), gr-qc/0003101. 

A. Dedner, F. Kemm, D. Kroner, C.-D. Munz, 
T. Schnitzer, and M. Wesenberg, J. Comput. Phys. 175, 
645 (2002). 

S. K. Godunov, Mat. Sb. 47, 271 (1959). 

M. J. Berger and J. Oliger, J. Comp. Phys. 53, 484 

(1984). 

L. Lehner, S. L. Liebling, and O. Reula, Class. Quant. 
Grav. 23, S421 (2006), gr-qc/0510111. 

B. Strand, J. Comput. Phys. 110, 47 (1994). 

G. Calabrese, L. Lehner, D. Neilsen, J. Pullin, O. Reula, 
O. Sarbach, and M. Tiglio, Class. Quantum Grav. 20, 
L245 (2003), gr-qc/0302072. 

A. Harten, P. D. Lax, and B. van Leer, SIAM Rev. 25, 
35 (1983). 

C. -W. Shu and S. Osher, J. Comput. Phys. 77, 439 
(1988). 

M. Cecere, L. Lehner, and O. Reula, Comput. Phys. 

Commun. 179, 545 (2008), 0801.2140. 

G. Toth, J. Comput. Phys. 161, 605 (2000). 

O. Rinne, L. Lindblom, and M. A. Scheel, Class. Quant. 

Grav. 24, 4053 (2007). 

M. Tiglio, L. Lehner, and D. Neilsen, Phys. Rev. D70, 
104018 (2004). 

G. Calabrese, L. Lehner, O. Reula, O. Sarbach, and 
M. Tiglio, Class. Quant. Grav. 21, 5735 (2004), gr- 
qc/0308007. 



13 



[49] L. Lehner and O. M. Moreschi, Phys. Rev. D76, 124040 

(2007), 0706.1319. 
[50] L. Baiotti, R. De Pietri, G. M. Manca, and L. Rezzolla, 

Phys. Rev. D75, 044023 (2007), astro-ph/0609473. 
[51] M. Berger and P. Collella, J. Comput. Phys. 82 82, 64 

(1989). 

[52] G. B. Cook, S. L. Shapiro, and S. A. Teukolsky, Astro- 
phys. J. 424, 823 (1994). 

[53] J. A. Font, T. Goodale, S. Iyer, M. Miller, L. Rezzolla, 
E. Seidel, N. Stergioulas, W.-M. Suen, and M. Tobias, 
Phys. Rev. D65, 084024 (2002). 

[54] T. W. Baumgarte and S. L. Shapiro, Astrophys. J. 585, 
930 (2003), astro-ph/0211339. 

[55] I. Hawke, E. Schnetter, L. Baiotti, and L. Rezzolla, Com- 
put. Phys. Commun. 169, 374 (2005). 

[56] M. D. Duez, Y. T. Liu, S. L. Shapiro, M. Shibata, and 
B. C. Stephens, Phys. Rev. Lett. 96, 031101 (2006), 



astro-ph/05 10653. 
[57] G. M. Manca, L. Baiotti, R. De Pietri, and L. Rezzolla, 

Class. Quant. Grav. 24, S171 (2007), 0705.1826. 
[58] M. W. Choptuik, Phys. Rev. Lett. 70, 9 (1993). 
[59] C. Gundlach and J. M. Martin-Garcia, Living Rev. Rel. 

10, 5 (2007), 0711.4620. 
[60] M. W. Choptuik, E. W. Hirschmann, and S. L. Liebling, 

Phys. Rev. D55, 6014 (1997), gr-qc/9701011. 
[61] S. L. Liebling, Phys. Rev. D58, 084015 (1998), gr- 

qc/9805043. 

[62] S. Noble, Ph.D. thesis, The University of Texas at Austin 
(2003). 

[63] K.-J. Jin and W.-M. Suen, Phys. Rev. Lett. 98, 131101 

(2007), gr-qc/0603094. 
[64] M. B. Wan, K. J. Jin, and W. M. Suen (2008), 0807.1710. 



